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INTRODUCTION 


The  feasibility  of  employing  heat  pipes  to  cool  the  hot  sections  of  the  Army's  ground-to- 
ground  missile  fins  has  been  studied.  A  portion  of  each  fin  is  exposed  directly  to  the  hot  gas  from  the 
exhaust  jet.  The  rest  of  the  fin  is  in  an  air  stream.  Heat  pipes  integrated  into  the  fin  structure  are 
proposed  for  transferring  heat  from  the  hot  section  to  the  cool  section  of  the  fin  to  prevent  damage 
to  the  fin  during  the  missile's  ten-minute  flight.  A  part  of  the  missile,  showing  the  location  of  fins  with 
respect  to  the  exhaust  jets,  is  depicted  in  Fig.  1.  The  fin  considered  in  this  study  is  a  thin  airfoil  and 
its  planform  is  rectangular.  The  fin  surface  is  symmetrical  with  respect  to  the  chord  line;  that  is,  the 
mean  chamber  line  is  coincident  with  the  chord  line.  The  angle  of  attack  is  zero  so  that  the  free- 
stream  flow  is  parallel  to  the  chord  line. 

The  proposed  fin  cooling  concept  consists  of  two  liquid-metal  heat  pipes  separated  by  a  planar 
rib  as  shown  in  Fig.  2.  Liquid  metal  has  been  considered  as  the  working  fluid  for  both  heat  pipes 
because  high  operating  temperatures  are  anticipated.  However,  the  use  of  liquid  metal  creates  a 
startup  problem  since  the  working  fluid  is  initially  in  a  solid  state  at  ambient  temperature.  These  heat 
pipes  thus  need  to  be  designed  so  that  they  start  successfully  from  the  frozen  state.  For  this  purpose, 
a  simple  model  based  on  the  lumped  heat-capacity  method  has  been  developed  for  predicting  the 
startup  and  transient  characteristics.  The  surface  of  the  fore  heat  pipe  in  the  leading  edge  region  has 
much  higher  heat  flux  than  the  heat  pipe  located  in  the  rear.  Therefore,  in  addition  to  the  transient 
analysis,  a  steady-state  analysis  has  been  performed  for  the  fore  heat  pipe  to  determine  whether  the 
heat  pipe  will  operate  as  designed. 

In  order  to  determine  the  heat  load  to  the  evaporator  and  the  heat  removal  rate  from  the 
condenser  for  the  fore  and  aft  heat  pipes,  it  is  necessary  to  evaluate  the  heat  transfer  coefficients 
between  the  surface  of  the  airfoil  fin  and  the  heat  source  and  sink.  A  fairly  good  approximation  may 
be  made  by  use  of  the  heat  transfer  correlations  for  smooth  cylinders  and  flat  plates  for  the  fore  and 
aft  heat  pipes,  respectively.  The  fore  heat  pipe  may  be  geometrically  transformed  to  a  circular 
cylinder  of  equivalent  diameter.  The  aft  heat  pipe  following  the  equivalent  cylinder  can  be  analyzed 
using  the  flat  plate  correlations. 

There  are  two  free-stream  flows:  exhaust  gas  stream  and  air  stream.  The  Mach  number  and 
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Fig.  1  Schematic  of  exhaust  jets  and  fins  of  the  missile. 
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static  temperature  of  each  stream  are  indicated  in  Fig.  3.  Since  the  exhaust  gas  flow  is  supersonic, 
a  curved  shock  wave  exists  in  front  of  the  leading  edge.  To  simplify  the  analysis,  a  normal  shock  is 
assumed  instead  of  the  curved  shock.  After  the  normal  shock,  the  Mach  number  and  static 
temperature  of  the  gas  become  M= 0.8892  and  r„e=1084K  as  can  be  calculated  from  the  normal 
shock  tables  (John,  1984). 

FORCED  CONVECTION  HEAT  TRANSFER  IN  HIGH-SPEED  FLOW 

The  startup  and  steady-state  characteristics  of  liquid-metal  heat  pipes  depend  on  the  heat 
source  and  heat  sink  conditions.  For  the  conditions  of  interest,  radiation  heat  transfer  has  been 
considered  and  found  negligible  compared  to  the  aerodynamic  heating.  Therefore,  only  aerodynamic 
heating  and  cooling  heat  transfer  rates  are  considered  here.  Considering  heat  transfer  on  the  surface 
of  heat  pipes  moving  at  high  velocities  in  the  hot  exhaust  gas  and  air  stream,  the  effects  of  the  fluid 
compressibility  and  viscous  dissipation  must  be  taken  into  account. 

The  heat  transfer  rate  Q  on  the  surface  in  high-speed  flow  can  be  calculated  with  the  same 
relations  used  for  low-speed  flow  if  the  average  heat  transfer  coefficient  h  is  defined  as 

Q  -  -  Tr)  (1) 

where  A  is  the  surface  area,  Tp  is  the  surface  temperature,  and  Tr  is  the  recovery  temperature.  The 
recovery  temperature  is  expressed  by  a  dimensionless  parameter,  r,  called  the  temperature  recovery 
factor  and  defined  as 


or 
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.  3  Flow  conditions  of  free  streams. 


r 


7  -  7 

r  “ 

7  -  7 

O  00 


(3) 


where  Tm  is  the  free-stream  temperature,  u  is  the  free-stream  velocity,  cp  is  the  specific  heat  at 
constant  pressure  of  the  fluid,  and  Ta  is  the  stagnation  temperature. 

Stagnation  temperature  for  a  perfect  gas  with  constant  specific  heats  is  found  from 


7  =  7  (1 

O 


u 


2c  7 

p  00 


(4) 


Expressing  this  equation  in  terms  of  the  Mach  number  gives 
T  -  T  (1  +  M 2) 

O  CD  V  OO  J 


(5) 


where  y  is  the  ratio  of  specific  heats  of  the  fluid.  Experiments  with  flat  plates  have  shown  that  the 
recovery  factor  r  can  be  approximately  expressed  in  terms  of  the  Prandtl  number  Pr.  For  laminar 
flow, 

r  =  Pr  Vl  (6) 

whereas  for  turbulent  flow, 

r  =  PrVa  ( 7 ) 

In  order  to  account  for  fluid  property  dependence  on  temperature,  Eckert  (1960)  has 
recommended  that  the  properties  be  evaluated  at  a  reference  temperature  T  given  by 

T'  =  Tm  ♦  0.5(7  -  7J  ♦  0.22(7  -  7J  (8) 

The  last  term  in  Eq.  (8)  can  also  be  expressed  by  the  Mach  number  as 

T’  ■  T.  .  0.5(7,  -  7.)  .  0.11(Y-l)rA/X  (9) 
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The  leading  edge  contour  may  be  approximated  by  a  cylindrical  shape  which  has  the  same 
circumference.  Over  the  forward  portion  of  the  cylinder  (0°<  Q  <80°)  where  the  angle  <9  is  measured 
from  the  stagnation  point,  the  empirical  equation  for  the  local  heat  transfer  coefficient  he  at  d  is  given 
as 


v> 

k 


1.14  RelD2Pr 


2/5 


3 

1- 

,90; 

(10) 


where  k  is  the  thermal  conductivity  of  fluid  and  ReD  is  the  Reynolds  number  based  on  the  diameter 
D.  That  is, 


(11) 


where  p  and  u  are  the  density  and  viscosity  of  the  fluid.  Equation  (10)  was  originally  presented  by 
Martinelli  et  al.  (1943)  for  the  range  of  6  from  0°  to  90°.  However,  Kreith  and  Bohn  (1986) 
suggested  to  use  this  equation  from  0°  to  80°,  as  has  been  done  in  the  present  work.  Integrating  Eq. 
(10)  along  the  arc  length  from  0  to  2nD/9  with  respect  to  D  0/2,  the  average  heat  transfer  coefficient 
h  becomes 


h  D 
k 


0. 94 Re p2Pr  25 


(12) 


The  heat  transfer  rate  on  the  surface  following  the  leading  edge  region  may  be  evaluated  by 
treating  the  surface  as  two  flat  plates.  The  flow  is  assumed  turbulent  over  the  entire  surface  because 
of  the  expected  free-stream  turbulence.  The  local  heat  transfer  coefficient  hx  for  the  turbulent 
boundary  layer  on  a  flat  plate  is  given  by  (Chapman,  1984) 


Nux- 


0.0296  RefPr  1/3 


(13) 


when  the  Reynolds  number  ranges  from  5x1 05  to  107  which  encompasses  the  range  of  flow  conditions 
considered  in  this  study.  Here,  the  characteristic  length  x  is  the  distance  starting  from  the  stagnation 
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point.  Integrating  Eq.  (13)  along  the  plate  length  L  from  xa  to  obtain  the  average  heat  transfer 
coefficient  h  gives 


Nul 


-  =  0.037  Re 
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l  \ 
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STEADY-STATE  ANALYSIS 


(14) 


The  leading  edge  region  of  the  missile  fin  is  subjected  to  a  severe  aerodynamic  heating 
environment.  Therefore,  the  fore  heat  pipe  must  be  able  to  handle  the  high  heat  flux  by  distributing 
the  heat  around  the  whole  circumference  of  the  wick  structure.  For  this  purpose,  it  is  desirable  to  use 
a  wick  structure  possessing  circumferentially  interconnecting  pores.  As  an  example,  the  present  study 
analyzes  a  screen  wick  saturated  with  liquid  sodium  for  characterizing  the  heat  pipe  performance  at 
steady  state.  The  operating  condition  of  the  fore  heat  pipe  must  be  below  the  operating  limits  on  its 
heat  transport  capability. 

Since  the  operating  temperature  and  the  condenser  length  of  the  heat  pipe  required  to  remove 
the  heat  load  are  not  known  a  priori,  iterative  calculations  are  necessitated.  First,  the  surface 
temperature  of  the  evaporator  Tp  e  is  assumed.  From  the  aerodynamic  heating  correlations,  the  heat 
load  to  the  evaporator  can  then  be  found.  Because  at  steady  state  the  total  heat  input  to  the 
evaporator  must  be  rejected  through  the  condenser  surface  to  the  heat  sink,  the  required  length  of  the 
condenser  can  be  determined.  Once  the  operating  temperature  of  the  heat  pipe  is  calculated,  all  the 
operating  limits  on  the  heat  transport  capability  are  evaluated  at  that  operating  temperature.  If  the 
operating  condition  exceeds  any  operating  limit,  a  new  value  of  the  surface  temperature  of  the 
evaporator  is  assumed  and  the  above  calculating  procedure  must  be  repeated  until  successful 
operation  is  obtained. 

The  heat  transfer  rate  to  the  cylindrical  surface  of  the  evaporator,  subtended  by  the  angle  6 
=  8k/9,  can  be  obtained  from 
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Q  =  h 
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(15) 


—  71  Dz 
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where  z£  is  the  given  evaporator  length  and  the  subscript  e  denotes  the  evaporator.  The  average  heat 
transfer  coefficient  for  the  evaporator  surface,  he ,  in  Eq.  (15)  can  be  calculated  from  Eq.  (12). 
Assuming  that  the  heat  applied  to  a  portion  of  the  circumference  of  the  heat  pipe  may  be  distributed 
equally  through  the  entire  heat  pipe  periphery  at  steady  state,  the  radial  heat  flux  can  be  expressed 
as 


= 


Qe 

%Dz 

C 


(16) 


This  assumption  can  be  made  because  liquid  metal  has  a  thermal  conductivity  large  enough  to  spread 
the  partially  applied  heat  around  the  circumference. 

It  is  important  to  determine  the  operating  temperature  of  the  heat  pipe.  The  temperature 
difference  between  the  vapor  and  the  liquid  at  the  liquid-vapor  interface  is  generally  small  and  can  be 
neglected.  Furthermore,  the  temperature  of  the  vapor  core  may  be  assumed  to  be  uniform  because 
of  the  small  thermal  resistance  for  vapor  flow.  Accordingly,  the  primary  temperature  drops  in  the 
heat  pipe  occur  through  the  pipe  wall  and  the  liquid-saturated  wick. 

Figure  4  shows  the  locations  where  temperatures  have  been  evaluated  and  the  dimensions  of 
the  heat  pipe.  Applying  the  Fourier  law  of  heat  conduction,  the  individual  temperature  drops  can  be 
written  as  follows  : 
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Fig.  4  Schematic  of  the 
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In  the  above  equations,  the  subscripts  e  and  c  indicate  the  evaporator  and  the  condenser  and  kp  is  the 
thermal  conductivity  of  pipe  material.  The  temperatures  Tp,  T^,  and  Tv  are  the  temperature  of  the 
pipe  wall  surface,  the  temperature  at  the  pipe-wick  interface,  and  the  vapor  temperature,  respectively. 
The  radii  ra ,  r, ,  and  rv  are  for  the  outside  and  inside  of  the  pipe  and  the  vapor  core.  The  effective 
thermal  conductivity  ke  for  the  liquid-saturated  screen  wick  in  Eqs.  (18)  and  (19)  can  be  found  from 

*z[(*z  ♦  K)  -  (1 -e)(*,  -  kj] 

k  =  _ _ i _ - : - —  (21) 

e  (kt  *  kj  *  (l-e)(kl  -  kj 


where  e  is  the  porosity  of  the  wick  structure,  k,  is  the  thermal  conductivity  of  the  liquid,  and  kw  is  the 
thermal  conductivity  of  the  wick  material. 

Notice  that  the  required  length  of  the  condenser  zc  in  Eqs.  (19)  and  (20)  to  remove  the  applied 
heat  load  is  unknown.  This  length  must  be  determined  such  that  the  operating  temperature  of  the 
heat  pipe  is  within  a  safe  range  to  avoid  damage  to  the  heat  pipe. 
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During  steady-state  operation,  all  the  heat  applied  to  the  evaporator  will  be  rejected  through 
the  surface  of  the  condenser.  The  equation  used  for  the  aerodynamic  heating  on  the  evaporator 
surface  can  also  be  employed  for  the  aerodynamic  cooling  on  the  condenser  surface  by  changing  the 
sign.  The  heat  removal  rate  from  the  condenser  can  then  be  found  from 
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In  Eq.  (22),  hc  and  Trc  require  the  reference  temperature  T  which  in  turn  requires  knowledge 
of  71  .  In  order  to  evaluate  Tpc ,  the  unknown  zc  in  Eqs.  (19)  and  (20)  must  be  assumed  so  that  Eq. 
(22)  is  satisfied.  If  the  heat  output  from  Eq.  (22)  is  not  equal  to  the  heat  input  from  Eq.  (15),  then 
a  new  value  of  zc  has  to  be  assumed  and  the  above  calculation  procedure  needs  to  be  repeated  until 
a  satisfactory  check  is  made. 

Using  the  described  technique,  the  operating  condition  of  the  heat  pipe  can  be  analyzed.  In 
order  to  confirm  whether  the  heat  pipe  operates  satisfactorily  under  this  condition,  the  maximum  heat 
transport  capability  of  the  heat  pipe  at  the  operating  temperature  must  be  evaluated  to  see  if  it 
exceeds  the  design  condition.  There  are  several  operating  limits  on  the  heat  pipe:  sonic,  entrainment, 
capillary,  and  boiling  limits.  Among  these  operating  limits,  the  lowest  one  at  a  given  operating 
temperature  provides  the  maximum  possible  value  of  heat  transfer  rate  at  that  temperature.  The 
boiling  limit  is  usually  not  a  problem  with  liquid-metal  heat  pipes  because  of  the  high  thermal 
conductivity  of  the  fluid  and  the  high  superheat  needed  to  initiate  boiling.  The  most  common  heat 
transfer  limit  is  the  capillary  limit.  The  capillary  limit  presented  here  is  for  the  most  extreme 
condition;  that  is,  the  heat  pipe  is  in  a  vertical  position  with  the  evaporator  above  the  condenser.  The 
sonic,  entrainment,  and  capillary  limits  can  be  calculated  from  the  expressions  below  (Chi,  1976). 
Sonic  limit: 

When  the  heat  pipe  is  operating  at  low  vapor  densities,  the  vapor  velocity  may  reach  the  speed 
of  sound.  For  this  situation  the  heat  transport  capability  of  the  heat  pipe  is  restricted  to  a  maximum 
value  by  the  existence  of  shocked  flow  at  the  evaporator  exit.  This  restriction  on  the  heat  transport 
capability  is  known  as  the  sonic  limit. 

Assuming  negligible  frictional  effects,  the  vapor  flow  analysis  gives  the  sonic  limit  Qs  as 
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where  A v  is  the  cross-sectional  area  of  the  vapor  core,  p„  is  the  vapor  density,  /i/g  is  the  latent  heat  of 
vaporization,  and  R  is  a  gas  constant. 

Entrainment  limit: 

During  the  heat  pipe  operation  the  liquid  and  the  vapor  flow  in  opposite  directions.  If  the 
vapor  velocity  is  sufficiently  high,  a  shear  force  exerted  by  the  vapor  at  the  liquid- vapor  interface  may 
entrain  droplets  of  liquid  and  carry  them  to  the  condenser  end.  When  this  happens,  the  wick  in  the 
evaporator  is  depleted  of  liquid  needed  for  continuous  evaporation.  This  limits  the  maximum  axial 
heat  transport  capability  and  is  called  the  entrainment  limit. 

The  entrainment  limit  Qem  is  determined  from 
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where  o  is  the  surface  tension  of  liquid,  and  rks  is  the  hydraulic  radius  of  the  wick  surface  pores.  The 
hydraulic  radius  of  the  screen  wick  can  be  found  from 


1 

rKs  -  ~ 


N 


-  d 


(25) 


where  N  and  d  are  the  mesh  number  and  the  wire  diameter  of  the  screen. 

Capillary  limit: 

The  maximum  capillary  pressure  must  be  greater  than  the  total  pressure  drop  associated  with 
the  liquid  and  vapor  flows  in  the  heat  pipe  for  a  normal  operation.  Otherwise,  the  wick  will  dry  out 
in  the  evaporator  due  to  insufficient  capillary  pumping  capability  to  return  the  liquid  to  the  evaporator 
from  the  condenser.  This  operating  limit  is  referred  to  as  the  capillary  limit. 

When  the  liquid  and  vapor  flows  are  laminar  and  incompressible  with  uniform  heat  fluxes 
applied  to  the  evaporator  and  condenser  sections,  the  capillary  limit  Qcap  can  be  found  from 
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In  the  above  equations,  fx,  and  fxv  are  the  liquid  and  vapor  viscosities,  K  is  the  permeability  of  the  wick 
structure,  and/v  and  Rev  are  the  Fanning  friction  factor  and  the  Reynolds  number  of  the  vapor  flow. 

Because  of  the  expected  high  operating  temperature,  a  sodium  heat  pipe  has  been  considered 
for  the  fore  heat  pipe.  Various  physical  properties  of  sodium  are  summarized  in  Table  I.  The 
container  material  is  type  304  stainless  steel  and  the  wick  is  four  layers  of  200  mesh  screen  of  type 
304  stainless  steel.  A  summary  of  the  sodium  heat  pipe  is  given  in  Table  II.  With  the  given  free- 
stream  conditions  shown  in  Fig.  3,  a  heat  transfer  rate  of  ge=204W  has  been  determined.  The 
average  radial  heat  flux  is  <?r=13  W/cm2.  The  calculated  results  are  listed  with  the  operating  limits 
in  Table  III.  As  can  be  noticed,  all  the  limits  are  much  higher  than  the  required  heat  transfer  rate. 
Therefore,  it  may  be  concluded  that  the  sodium  heat  pipe  with  the  total  length  of  0.078  m  operates 
successfully  under  given  heat  source  and  sink  conditions. 
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Table  I.  Physical  Properties  of  Sodium 


Molecular  Weight 

22.991  - 

Melting  Point 

371. OK 

Boiling  Point  at  1  atm 

1151.2K 

Critical  Temperature 

2500K 

Critical  Pressure 

370  bar 

Critical  Density 

180  kg/m3 

Latent  Heat  of  Fusion 

113.044  kJ/kg 

Latent  Heat  of  Vaporization  at  1151. 2K 

3927.1  kJ/kg 

Liquid  Thermal  Conductivity  at  900K 

61.4  W/m-K 

Vapor  Thermal  Conductivity  at  900K 

4.06  X  10’2  W/m-K 

Liquid  Viscosity  at  900K 

2.02  X  10-4  N-s/m2 

Vapor  Viscosity  at  900K 

2.06  X  10'5  N-s/m2 

Surface  Tension  at  900K 

0.142  N/m 

Table  II.  Fore  Heat  Pipe  Details 

Working  Fluid 

Sodium  (Na) 

Container  Material 

Type  304  Stainless  Steel 

Total  Heat  Pipe  Length 

0.078  m 

Evaporator  Length 

0.05  m 

Condenser  Length 

0.028  m 

Container  Outside  Diameter 

0.01m 

Container  Inside  Diameter 

0.007  m 

Wick  Material 

Type  304  Stainless  Steel 

Wick  Structure 

Four  Layers  of  200  Mesh  Screen 
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Table  III.  Calculated  Results  for  the  Fore  Heat  Pipe 


Qe  =  204W 

qr  =  13  W/cm 

re  =  0.8526 

rc  =  0.8289 

Tre  =  1211.7K 

Tr>c  =  326.6K 

T;  =  1055. IK 

T*  =  625.5K 

Tp>e  =  970K 

Tpc  =  939.3K 

Tpw,e  -  960.7K 

Tmc  =  956K 

Tv  =  959K 

Sonic  Limit,  Qs  =  1764.9  W 
Entrainment  Limit,  Qen,  =  1121.3W 
Capillary  Limit,  Qcap  =  441  W 


It  is  of  importance  to  notice  that  the  steady-state  analysis  mentioned  above  is  restricted  to  the 
situation  when  there  is  no  free  molecular  region  in  the  heat  pipe.  If  the  heat  pipe  is  steadily  operating 
with  a  portion  of  the  free  molecular  region  in  the  inactive  condenser  section,  the  results  from  the 
steady-state  analysis  will  be  in  error.  As  explained  in  the  next  section,  the  transient  model  predicts 
the  existence  of  the  free  molecular  region. 

TRANSIENT  ANALYSIS 

When  the  liquid-metal  heat  pipe  is  started  up  from  ambient  temperature,  the  working 
substance  is  in  the  solid  state.  The  vapor  density  is  so  small  that  free  molecular  flow  regime  prevails 
over  the  entire  vapor  core.  As  heat  is  applied  to  the  evaporator,  the  frozen  working  substance  will 
be  melted.  As  heating  continues,  the  free  molecular  flow  in  the  evaporator  will  become  continuum 
flow.  Evaporation  can,  thereafter,  take  place  at  the  liquid-vapor  interface  resulting  in  the  increase 
of  the  vapor  temperature  of  the  evaporator.  With  further  heating,  a  continuum  flow  front  moves 
toward  the  condenser  end  while  the  sonic  limit  exists  at  the  evaporator  exit.  Once  the  front  has 
reached  the  end  of  the  condenser,  continuum  flow  exists  over  the  entire  heat  pipe  length.  The 
operating  temperature  will  then  continue  to  rise  until  steady  state  is  reached.  The  solution  procedure 
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is  the  same  for  the  startup  prediction  of  the  fore  heat  pipe  and  the  aft  heat  pipe  except  for  the 
calculations  of  the  heat  transfer  coefficients,  he  and  hc ,  and  the  heat  input  and  output. 

The  transition  from  the  free  molecular  flow  regime  to  the  continuum  flow  regime  occurs  when 
the  Knudsen  number  Kn  is  less  than  0.01. 


Kn 


=  —  <  0.01 
/ 


(30) 


Here,  X  is  the  mean  free  path  of  vapor  and  /  is  a  characteristic  length.  From  kinetic  theory,  the 
dynamic  viscosity  /uv  and  the  mean  molecular  velocity  v  are 

pv=  0.499  pyvA  (31) 


and 
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From  Eqs.  (30)  -  (32),  it  can  be  noticed  that  continuum  flow  exists  when  the  temperature  of  the  vapor 
reaches  the  transition  temperature  T,  given  by 
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Iterations  are  required  to  determine  T,  because  of  the  temperature  dependence  of  piv  and  pv.  The 
characteristic  length  /  in  Eq.  (33)  is  the  vapor  core  diameter  for  the  fore  heat  pipe  and  is  assumed  to 
be  the  maximum  thickness  of  the  vapor  core  for  the  aft  heat  pipe. 

The  lumped  heat-capacity  method  has  been  applied  to  the  heat  pipe.  For  the  startup  from  the 
frozen  state,  the  evaporator  temperature  can  be  found  from  the  energy  balance  on  the  evaporator 
section  as 
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(34) 


when  Tp<Tt .  Here,  the  superscripts  n  and  n+1  denote  times  t  =  nAt  and  t  =  (n+l)At  where  At  is 
a  time  increment.  The  effective  volumetric  heat  capacity  per  unit  spanwise  length  C  depends  on 
whether  the  lumped  temperature  T”  is  greater  or  less  than  the  melting  temperature  Tm  of  the  working 
substance.  When  T”<Tm, 

c  ■  (Pcp)FJr  *  •  e(p cp)A„  (35) 


and  when  Tp">Tm, 


C  "  ^Cp)pAp  +  (l-eKPcp)wA> 
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The  subscripts/?,  w,  s,  and  l  denote  the  pipe  wall  material,  wick  material,  and  the  solid  and  liquid  state 
of  the  working  substance.  Ap  and  are  the  cross-sectional  areas  of  the  pipe  wall  and  the  wick 
structure.  In  Eqs.  (35)  and  (36)  the  heat  capacity  of  the  vapor  core  is  not  included  since  it  is 
negligible. 

When  T”  is  greater  than  the  transition  temperature  T, ,  the  lumped  temperature  Tpn+l  after  a 
time  increment  At  can  be  evaluated  depending  on  the  location  of  the  continuum  front  z": 
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Q{  in  Eq.  (38)  is  the  heat  transferred  from  the  continuum  region  to  the  control  volume  as  indicated 
in  Fig.  5  and  is  given  by 


2,-2, 
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Here,  hs,  is  the  latent  heat  of  fusion  of  the  working  substance. 

The  sonic  limit  Qs  in  Eq.  (40)  can  be  evaluated  from  Eq.  (23).  Because  of  the  complex  cross- 
sectional  shape  of  the  fin,  the  vapor  core  area  Av  of  the  aft  heat  pipe  for  Eq.  (23)  has  been  evaluated 
by  Simpson's  rule.  For  the  fore  heat  pipe,  the  heat  transfer  rates  are  calculated  from  Eqs.  (15)  and 
(22)  with  the  heat  transfer  coefficient  from  Eq.  (12).  For  the  aft  heat  pipe,  the  heat  transfer 
coefficient  for  the  portion  of  the  fin  following  the  equivalent  cylinder  is  calculated  from  the  flat  plate 
equation.  Using  the  average  heat  transfer  coefficient  from  Eq.  (14),  the  heat  transfer  rate  to  the 
heated  and  the  cooled  zones  of  the  aft  heat  pipe  can  be  calculated  from 


Qe  -  2 he(L-xo)ze(Trl 
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and 


2 hc(L-  xo)(zn 


(42) 


The  location  of  the  continuum  front  can  be  found  by  applying  the  energy  balance  to  the  small 
control  volume  of  width  (z”+1  - 1)  as  shown  in  Fig.  5. 
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where  the  effective  volumetric  heat  capacity  Cf  of  the  free  molecular  flow  region  depends  on  whether 
the  temperature  of  the  free  molecular  flow  region  Tfc  is  greater  or  less  than  the  melting  temperature 
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Fig.  5  Propagation  of  the  continuum  front. 


Tm  as  in  Eqs.  (35)  and  (36). 

A  computer  program  has  been  developed  for  predicting  transient  startup  of  liquid-metal  heat 
pipes.  A  listing  of  a  Fortran  version  of  this  program  is  given  in  Appendix.  The  Fortran  variable 
names  assigned  to  the  various  quantities  are  presented  in  the  first  part  of  the  program.  The  computer 
program  consists  of  a  main  program  and  three  subroutines.  These  subroutines  are  for  the 
thermophysical  properties  of  air,  type  304  stainless  steel,  and  sodium.  The  properties  are  excerpted 
from  Brennan  and  Kroliczek  (1979),  Incropera  and  DeWitt  (1981),  and  Vargaftik  (1975). 

The  transient  model  has  been  applied  to  the  fore  and  aft  heat  pipes  to  investigate  the  feasibility 
of  using  the  liquid-metal  heat  pipe.  Sodium  is  the  working  fluid  and  type  304  stainless  steel  is  used 
for  the  container  and  wick  structure.  Descriptions  of  the  fore  and  aft  heat  pipes  are  given  in  Tables 
II  and  IV. 

The  computed  results  for  the  fore  heat  pipe  are  plotted  in  Figs.  6  and  7.  These  figures  show 
the  variation  of  the  lumped  temperature  and  the  continuum  flow  front  as  a  function  of  time.  Once 
the  temperature  exceeds  the  transition  temperature  of  725. 2K,  the  continuum  front  starts  to 


Table  IV.  Input  Data  for  Startup  Prediction  for  the  Aft  Heat  Pipe 


Working  Fluid 
Container  Material 
Plate  Length 

Total  Available  Heat  Pipe  Length 

Evaporator  Length 

Wall  Thickness 

Wick  Thickness 

Wick  Porosity 

Characteristic  Length 

Vapor  Core  Area 


Sodium  (Na) 

Type  304  Stainless  Steel 

0.12m 

0.25m 

0.05m 

0.0015m 

0.0005m 

0.7 

0.014m 

0.00067m2 
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propagate  from  the  evaporator  exit  toward  the  condenser  end.  When  the  total  spanwise  length  of  the 
fin  is  used  as  the  total  length  of  the  heat  pipe  (ze  +  zc  -  0.25  m),  steady  state  is  reached  at  65  seconds 
with  Tp  =  819K  and  z  =  0.08  m.  The  continuum  front  does  not  move  to  the  condenser  end  because 
the  operating  condition  of  the  heat  pipe  is  lower  than  the  design  condition.  A  shorter  length  than  z 
=  0.08  m  may  be  used  for  the  given  conditions  with  a  much  higher  operating  temperature;  then  the 
continuum  front  will  reach  the  condenser  end  and  a  normal  operation  can  be  obtained. 

Figures  6  and  7  include  the  computed  results  for  ze  +  zc  =  0.078  m  which  is  the  same  length 
as  that  in  the  steady-state  analysis.  The  front  reaches  the  condenser  end  at  40  seconds.  Thereafter, 
the  temperature  increases  until  steady  state  is  reached  at  120  seconds  with  Tp  =  985K.  This  steady 
operating  temperature  is  close  to  Tv  =  959K  from  the  steady-state  analysis. 

Similar  calculations  have  been  conducted  for  the  aft  heat  pipe.  The  transition  temperature  is 
689. IK,  and  the  results  are  shown  in  Figs.  8  and  9.  When  the  total  length  of  the  fin  is  used  as  the 
total  heat  pipe  length,  i.e.,  ze  +  zc  =  0.25  m,  the  continuum  front  stops  at  z  =  0.1115  m  with  Tp  = 
810K  after  120  seconds.  When  a  shorter  length  of  ze  +  zc  =  0.078  m  is  used,  the  continuum  front 
reaches  the  condenser  end  at  18  seconds.  Since  then  the  temperature  increases  until  steady  state  is 
reached  at  120  seconds  with  Tp  =  997K. 

With  a  heat  pipe  length  of  0.078  m,  the  operating  temperature  of  the  fore  heat  pipe  ( Tp  = 
985K)  is  slightly  lower  than  that  of  the  aft  heat  pipe  ( Tp  =  997K).  Even  though  very  high  heat  flux 
is  anticipated  at  the  fore  heat  pipe,  its  operating  temperature  is  lower  due  to  the  distributed  heat  input 
around  the  circumference  of  the  wick  structure. 

CONCLUSIONS  AND  RECOMMENDATIONS 

The  use  of  heat  pipes  to  cool  the  hot  region  of  a  missile  fin  has  been  studied.  A  steady -state 
analysis  has  been  performed  for  the  heat  pipe  in  order  to  determine  whether  the  operating  condition 
satisfies  the  limitations  on  the  heat  transport  capability  when  the  heat  pipe  is  normally  operating 
without  any  free  molecular  region.  In  addition,  startup  characteristics  of  the  heat  pipe  from  the  frozen 
state  have  been  analyzed  by  a  lumped  heat-capacity  method.  This  transient  model  provides  the 
necessary  length  of  the  heat  pipe  and  predicts  the  temperature  and  continuum  front  location  until 
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steady  state  is  reached.  Using  the  transient  analysis,  it  has  been  found  that  the  proper  length  of  the 
condenser  is  critical  for  heat  pipe  design.  If  the  condenser  length  is  too  short,  the  pipe  temperature 
increases  dramatically.  Meanwhile,  if  it  is  too  long,  there  is  a  risk  of  freezing  the  vaporized  sodium 
in  the  inactive  condenser  region,  causing  dryout  in  the  evaporator  due  to  insufficient  liquid. 

Since  the  current  transient  model  uses  the  lumped  heat-capacity  method,  it  cannot  predict 
whether  dryout  occurs  in  the  evaporator  due  to  an  insufficient  amount  of  liquid  to  replenish  the 
evaporator  because  of  a  suddenly  applied  heat  load.  Since  a  non-condensable  gas  helps  the  startup 
of  the  liquid-metal  heat  pipe  from  the  frozen  state,  it  is  often  desirable  to  include  a  little  non¬ 
condensable  gas  in  the  vapor  space.  However,  in  the  present  case  it  is  not  possible  to  use  non¬ 
condensable  gas  if  the  missile  is  severely  maneuvering  so  that  the  orientation  of  the  heat  pipe  is 
changing.  The  reason  for  this  is  as  follows.  When  the  evaporator  is  located  below  the  condenser, 
the  non-condensable  gas  must  have  a  smaller  molecular  weight  than  that  of  the  working  fluid. 
Meanwhile,  when  the  evaporator  is  above  the  condenser,  the  molecular  weight  of  the  non¬ 
condensable  gas  must  be  larger  than  that  of  the  working  fluid.  Otherwise,  the  lighter  gas  will  be 
swept  to  the  evaporator  and  cover  the  wick  surface  causing  a  hot  spot.  The  hot  zone  may  then 
exceed  the  melting  point  of  pipe  material. 

It  is  suggested  that  the  wick  structure  for  the  fore  heat  pipe  consist  of  interconnecting  pores 
so  as  to  spread  the  localized  heat  to  the  peripheral  wick.  For  the  aft  heat  pipe,  however,  the  wick 
structure  with  interconnecting  pores  is  not  promising.  This  is  because  missile  maneuvering  may  cause 
the  melted  liquid  to  collect  in  the  trailing  edge  region.  The  most  reliable  wick  structure  may  be  axial 
grooves  implanted  with  sintered  metal  fibers.  The  metal  fibers  must  provide  a  sufficient  wicking 
height  for  heat  pipe  operation  against  gravity.  Each  groove  needs  to  be  partially  filled  with  metal 
fibers  so  that  there  is  no  cross  communication  with  liquid.  This  is  a  non-standard  wick  design  and 
involves  some  technical  risk. 

Both  the  steady-state  analysis  and  the  transient  model  rely  on  many  assumptions.  Also,  the 
free-stream  conditions  cannot  possibly  be  the  same  for  each  missile  flight.  Based  on  this,  it  is 
recommended  that  the  heat  pipe  be  made  with  a  slightly  extra  length  for  the  condenser  than  necessary 
so  that  some  free  molecular  region  is  allowed.  However,  in  order  to  confirm  the  feasibility  of  using 
the  liquid-metal  heat  pipe  to  cool  the  hot  region  of  the  missile  fin,  experiments  must  be  performed 
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under  both  transient  and  steady-state  conditions.  The  heat  pipe  also  needs  to  be  tested  at  various  tilt 
angles.  The  most  critical  tilt  condition  is  that  of  the  heat  pipe  in  a  vertical  position  with  the 
evaporator  above  the  condenser. 
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APPENDIX 

LISTING  OF  THE  COMPUTER  PROGRAM 
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program  HPWING 

common/blkl/ao (8) , al (8) , a2 (8) , a3 (8) , a4 (8) , a5  (8) 
common/blk2/airk,  airrhov, aircp, airmuv, airgam, airPr, airR 
common/blk3/ subrhov, subrhol, subrhos, submuv, subhfg, subgam, 

+  subhsl, subR, subcpl, subcps , subcpv, submul, subsig, subk 

common/blk4/piperho, pipecp, wickrho, wickcp 
common/blk5 / Tprop 
common/blk6/iwf 

character *20  filechk, filedat , filein, filepot, fileout, junk 

real  Len,Me,Mc 

real  Tp (1500) , z (1500) 

**  *************  ****************  ****************************************** 


c 

VARIABLE 

NAJ 

c 

c 

ao 

= 

c 

al 

= 

c 

a2 

a 

c 

a3 

= 

c 

a4 

= 

c 

a5 

a 

c 

c 

aircp 

= 

c 

airgam 

= 

c 

airk 

= 

c 

airmuv 

= 

c 

airPr 

= 

c 

airR 

a 

c 

airrhov 

= 

c 

Ap 

= 

c 

Av 

= 

c 

Aw 

= 

c 

cval 

= 

c 

cvalf 

a 

c 

delt 

a 

c 

diam 

a 

c 

diamlam 

a 

c 

epscp 

= 

c 

epsPr 

= 

c 

c 

filechk 

a 

c 

filedat 

= 

c 

filein 

= 

c 

filepot 

ss 

c 

fileout 

= 

c 

c 

hbarc 

a 

c 

hbare 

= 

c 

c 

istep 

= 

c 

c 

Len 

= 

c 

c 

Me 

= 

c 

Me 

a 

c 

Pi 

= 

c 

pipecp 

= 

c 

piperho 

= 

c 

por 

= 

c 

coefficient  in 
coefficient  in 
coefficient  in 
coefficient  in 
coefficient  in 
coefficient  in 


curve  fit 
curve  fit 
curve  fit 
curve  fit 
curve  fit 
curve  fit 


for  working  fluid 
for  working  fluid 
for  working  fluid 
for  working  fluid 
for  working  fluid 
for  working  fluid 


specific  heat  of  air  (J/kg  K) 
ratio  of  specific  heats  of  air 
conductivity  of  air  (W/m  K) 
viscosity  of  air  (N  s/mA2) 

Prandtl  #  of  air 
gas  constant  of  air  (kJ/kg-K) 
density  of  air  (kg/mA3) 
cross-sectional  pipe  area  (mA2) 

calculated  vapor  core  area  by  Simpson's  rule  (mA2) 
cross-sectional  wick  area  (mA2) 

effective  volumetric  heat  capacity  per  unit  length 
same  as  cval  except  for  the  free  molecular  region 
time  increment  (sec) 

equivalent  vapor  core  diameter  of  heat  pipe 
vapor  characteristic  length  used  for  Ttr  calculation 
accuracy-check  value  for  cp 
accuracy-check  value  for  Pr  # 

checklist  filename 

output  property  data  for  review 

design  input  data 

working  fluid  coefficient  filename 
output  filename 

avg  heat  transfer  coefficient,  cond 
avg  heat  transfer  coefficient,  evap 


iteration  step  # 

length  from  stagnation  point 
to  trailing  edge  along  surface 
Mach  #  of  flow  over  condenser 
Mach  #  of  flow  over  evaporator 
3.14159 

specific  heat  of  pipe  (J/kg-K) 
density  of  pipe 
porosity  of  wick 


(m) 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c* 


Qc 

Qe 

Qf 

Q1 

Qs 

re 

rc 

subcpl 

subcps 

subcpv 

sub gam 

subhfg 

subk 

submul 

submuv 

subR 

subrhol 

subrhos 

subrhov 

subsig 

Tf 

thkwall 

thkwick 

Tinfe 

Tinfc 

tmax 

Tmelt 

Tp 

Tprop 

Trecc 

Trece 

Trval 

Tstagc 

Tstage 

Tstarc 

Tstare 

Ttr 

wickcp 

wickrho 

xo 


z 

zc 

ze 

zt 


heat  removal  rate  from  condenser 
heat  load  to  evaporator 

heat  from  continuum  region  to  molecular  flow  region 
latent  heat  (melting  process) 
sonic  limit 

recovery  factor  in  evaporator  section 
recovery  factor  in  condenser  section 

liquid  cp  of  working  substance 

solid  cp  of  working  substance  (J/kg-K) 

vapor  cp  of  working  substance  (J/kg-K) 

ratio  of  specific  heats  of  working  substance 

heat  of  vaporization  of  working  substance 

conductivity  of  working  substance 

liquid  viscosity  of  working  substance 

vapor  viscosity  of  working  substance 

gas  constant  of  working  substance 

liquid  density  of  working  substance 

solid  density  of  working  substance 

vapor  density  of  working  substance 

surface  tension  of  working  substance 

temperature  in  free  molecular  region  (K) 
thickness  of  wall  (m) 
thickness  of  wick  (m) 

free  stream  temp  over  evaporator  (K) 
free  stream  temp  over  condenser  (K) 

max  time  to  be  considered  (sec) 
melt  temp  (K) 
pipe  temp  (K) 

temp  to  evaluate  properties  at  (K) 
recovery  temp,  condenser*  (K) 
recovery  temp,  evaporator  (K) 
increment  to  Ttr  to  iteratively  find  Ttr 
stagnation  temp,  condenser  (K) 
stagnation  temp,  evaporator  (K) 
reference  temp,  condenser  (K) 
reference  temp,  evaporator  (K) 
transition  temp  (K) 

specific  heat  of  wick  (J/kg-K) 
density  of  wick 

location  beyond  which  the  flat  plate 
correlation  applies  (m) 

position  along  flat  plate  (m) 

length  of  condenser  (m) 

length  of  flat  plate  (m) 

total  available  length  of  heat  pipe  (m) 


write (*, *) 
write (* , *) 
write (*, *) 
write (* , *) 


/  *  *' 

'*  Enter  Working  Fluid  *' 
'  *  *' 
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write (*,*)  ' *  1  =  Sodium  (Na)  *' 

write (*,*)  '*  2  =  Potassium  (K)  *' 

write (*,*)  '*  *' 

write  (*, *)  ' ************************** ' 
write (*, * ) 

read(*,*)  iwf 

if (iwf .eq. 1)  filepot-' sodium.dat' 

if ( iwf . eq . 2 )  filepot-' pot assum . dat ' 

fileout-' missile .out' 

f ilechk-' missile . chk' 

filedat-' missile .prp' 

filein  =' missile . inp' 

open (8/ file-filepot, status-' unknown' ) 

open (9, file-filedat , status-' unknown' ) 

open (10, f ile=f ilechk, status-' unknown' ) 

open (11 , file-fileout , status-' unknown' ) 

open (12, file-filein, status-' unknown' ) 

c*************************************************** 


pi-2 . *acos  (0) 

c  ***  Time  step  in  sec  *** 
c 

write (*,*)  'Enter  value  for  time  increment  :  (>0.5)' 

read (*, *)  delt 

c  ***  Read  in  Liquid/Vapor  Working  Fluid  Equation  Coefficients 

do  1  i-1,8 

read  (8,  *)  ao(i)  ,  al  (i)  ,  a2  (i)  ,  a3  (i)  ,  a4  (i)  ,  a5  (i) 

1  continue 


c 

c 

c 


Design  Input  Variables  *** 


Data  to 
read (12, 
read (12, 
read (12, 
read (12, 
read (12, 
read (12, 
read ( 12 , 
read (12, 
read (12, 
read (12, 
read (12, 
read (12, 
read(12, 
read (12, 
read(12, 
read (12, 
read (12 , 
read (12, 
read (12, 
read(12, 
write ( *, 
write (*, 
write (*, 
write (*, 


be  read  in 
1000, end-2) 

*)  Len 
1000, end-2) 

*)  ze 

1000, end-2) 

*)  zt 

1000, end-2) 

*)  diam 
1000, end-2) 

*)  diamlam 
1000, end-2) 

*)  thkwall 
1000, end-2) 

*)  thkwick 
1000, end-2) 

*)  por 
1000, end-2) 

*)  Av 

1000, end-2) 
tmax 

'All  input  DATA  has  been  read  in. 


from  file  missile. inp 
junk 

junk 

junk 

junk 

junk 

junk 

junk 

junk 

junk 

junk 


'Len 
'  ze 


-  ' , Len 

-  '  ,ze 


★  ★  ft  ★ 
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write (*, *) 
write (*, *) 
write  (*, *) 
write  (*, *) 
write (*, *) 
write (*, *) 
write  (*, *) 
write  (*, *) 
write (*, *) 
write  (11, *) 
write (11, *) 
write (ii, *) 
write (li, *) 

write (11, *) 
write (11, *) 
write (11, *) 
write (li,  *) 
write (li, *) 
write (11, *) 
write (li, *) 
goto  3 


'zt 

' diam  = 
'  diamlam= 
'  thkwall= 
'  thkwick= 
'por  = 
Av 

tmax  = 
'Len 

'ze  = 

'zt 

'  diam  = 
'  diamlam= 
f thkwall= 
r thkwick= 
por  = 
Av 

tmax  = 


'f  zt 
/  diam 
'  f diamlam 
/ thkwall 
r thkwick 
/Por 
/Av 
/  tmax 

'/Len 
'  /  ze 
['Zt 
/  diam 
/  diamlam 
, thkwall 
/ thkwick 
/P  or 
/Av 
/  tmax 


c 

2 


Lengths  are  in  meters 


Len  = 

ze  = 

Zt  = 

diam  = 

diamlam* 
thkwall= 
thkwick= 
por  = 
Av  - 

tmax  = 


:  0.12 
0.05 
0.25 
0.0622 
0.014 
0.002 
0 . 004 
0.7 

0 . 00067 
20. 


r3=diam/2 . 


Ap=pi*ri**2 . 
Aw=pi *r2  * *2 . 


tmax=tmax*60 


pi*r2**2 . 
pi*r3**2 . 


Working  Substance 

for  sodium: 
if(iwf.eq.i)  then 
Tmelt  =  370.85 
subR  »  8315. /23. 
endif 

for ^ potassium 
if (iwf ,eq.2)  then 
Tmelt  =  336.4 
subR  =  8315. /39.1 

UendJf°f  ab0V®  are  J/k94-K 


Flow  Conditions 


Tinfe  =  1084. 

Tinfc  »  299.8 

Me  =  0.8892 

Me  =0.75 

airR  =  287. 

c  units  of  above  are  J/kg-K 
game  =  1.4 
game  =  1.4 

c  Accuracy-Check  Values 

epsep  =  10. 
epsPr  =  0.01 
epsTr  =  5.00 

c  Set  temp  to  Tinfe  to  evaluate  properties  for  first  iteration 

c  Velocity  of  flow  (m/s) 

uc=Mc* (gamc*airR*Tinfc) **0 . 5 
ue=Me* (game* airR* Tinfe) **0.5 

c******  Begin  Calculations  ******** 
istep=l 
time*0 . 
z (istep) =0 . 

Tp (istep) =Tinfc 
c 

c  Want  to  consider  Ttr  to  be  a  f (T) /  therefore,  iterative 

c  Guess  Ttrnew  =  600.  as  a  start  point 


Ttrnew=600 . 
tprop=Ttrnew 
call  subprop 

50  Ttr=4  990  .  *pi/subR*  (submuv/  (subrhov*diamlam)  )  **2  . 

Trval=l . 0 

if (abs (Ttr-Ttrnew) .It. 100.)  Trval=0.01 

if (abs (Ttr-Ttrnew) .It. epsTr)  goto  90 

if (Ttr . It .Ttrnew)  then 

Ttrnew=Ttrnew-Trval 

tprop=Ttrnew 

call  subprop 

goto  50 

else 

Ttrnew*Ttrnew+Trval 
tprop=Ttrnew 
call  subprop 
goto  50 
endif 

90  write (*,*)  'Ttransition  =  ' , Ttr 

write  (11,*)  'Ttransition  -  '  ,  Ttr 
write (11, *) 

write  (11,*)  'Step  t  [s]  Temp  [K]  Front  [m]  he  [W/m/"2-K] 
+hc  [W/m/v2-K]  ' 

write  (11,*)  ' -  -  -  -  - 

A _ f 
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write  (11, 1030)  istep, time, tp (istep) , z (istep) 


Tprop*Tinfe 
call  airprop 

101  re=airPr** (1 . /3 . ) 

Tstage-Tinfe* (1 .+  ( (game-1 .) /2 .) *Me**2 . ) 

Trece=Tinfe+re* (Tstage-Tinfe) 

Tstare= (Tp (istep) +Tinfe) /2 .+0 . 22*re* (game-1 - ) /2 . *Me**2 . *Tinfe 

c 

c  Now  evaluate  cp  and  Pr  at  Tstare;  compare  to  assumed  value, 

c 

cpold=aircp 

Prold=airPr 

Tprop^Tstare 

call  airprop 

cpnew^aircp 

Prnew=airPr 

dif fcp=abs (cpnew-cpold) 
dif fPr=abs (Prnew-Prold) 

if  (diffcp.gt .epscp  .or.  diffPr .gt .epsPr)  goto  101 

c  If  OK,  then  evaluate  hbar,  else  with  new  Pr  go  back  and  calc  r 

c  (actually  go  back  to  101  and  calculate  based  on  properties  at  T*) 

c  write  (10,*)  'Pr#  and  cp  converged  to  '  , Prnew, cpnew 

c  write  (10, 1010)  istep, Tp (istep) , Tstage, Trece, Tstare 

c  write  (10, *)  ' airmuv,  Len  =  ',airmuv,Len 

c  write (10,*)  'airk,  airrhov  »  9  , airk, airrhov 

c  write (10,*)  'ue,  airPr  «  ',ue,airPr 

xo=pi*0 . 01*2 . /9 . 
xl=xo/Len 

hbare=airk/Len*0 . 037* (airrhov*ue*Len/airmuv) **0 . 8*airPr** ( 1 . / 3 . ) 
+  * (1 . -xl**0 . 8) /  (l.-xl) 

c  write (10,*)  ' hbare  =  ',hbare 

c  --  Now  do  same  process  for  condenser  - 

Tprop^Tinfc 
call  airprop 

107  rc=airPr** (1 . /3 . ) 

Tstagc=Tinfc*  (1 .+  (  (gamc-1 . ) /2 . )  *Mc**2 . ) 

Trecc=Tinfc+rc* (Tstagc-Tinfc) 

Tstare25  (Tp  (istep)  +Tinfc)  / 2  .  +0 .22*rc*  (gamc-1 . )  /2  .  *Mc**2  .  *Tinfc 
c 

c  Now  evaluate  cp  and  Pr  at  Tstare;  compare  to  assumed  value, 

c 

cpold=aircp 

Prold=airPr 

Tprop^Tstarc 

call  airprop 

epnew^airep 

Prnew«airPr 

diffcp=abs (cpnew-cpold) 
dif fPr=abs (Prnew-Prold) 

if  (diffcp.gt .epscp  .or.  dif fPr .gt .epsPr)  goto  107 
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If  OK,  then  evaluate  hbar,  else  with  new  Pr  go  back  and  calc  r 
(actually  go  back  to  107  and  calculate  based  on  properties  at  T*) 

write(10,*)  'Pr#  and  cp  converged  to  ' , Prnew, cpnew 
write  ( 10 , 1011)  istep, Tp (istep) , Tstagc, Trecc, Tstarc 

hbarc=airk/Len*0 .037* (airrhov*uc*Len/ainauv) **0 . 8*airPr**  (1 . /3 . ) 
*  (1 . -xl**0 . 8) / (l.-xl) 
write  (10,*)  'hbarc  =  ',hbarc 

Tprop-Tp (istep) 
call  subprop 
call  pipprop 

Qs »subrhov*subhfg*Av*  ( (subgam*subR*Tp  (istep)  )  / 

(2 . * (subgam+1 . ) ) ) **0.5 
Qe=2 . *hbare* (Len-xo) *ze* (Trece-Tp (istep) ) 

Qc=2 . *hbarc* (Len-xo) * (z (istep) -ze) * (Tp (istep) -Trecc) 
if (Qc . It . 0 . )  Qc=0 . 

write  (10, 1020)  Qs, Qe, Qc, subrhov, subhfg 
if (Tp (istep) .gt.Ttr)  then 

cval=*  (piperho*pipecp*Ap)  +  (1 .  -por)  *  (wickrho*wickcp*Aw)  + 
por* (subrhol*subcpl*Aw) 

if (z (istep)  . ge. ze.and.z (istep)  .It .  (zt) )  then 

Tp (istep+1) ^Tp (istep)  +  (Qe-Qc-Qf ) *delt/ (cval*z (istep) ) 
Qf=Qs-Qc- (cval* (z (istep) -ze) * (Tp (istep+1) -Tp (istep) ) /delt) 
Quse^subrhol *por *Aw* subhs 1 
Qf=Qf-Quse* (z (istep) -z (istep-1) ) 

Tprop=Tinfc 
call  subprop 
call  pipprop 

cvalinf= (piperho*pipecp*Ap)  +  (1 . -por) * (wickrho*wickcp*Aw)  + 
por* (subrhos*subcps*Aw) 

Tprop^Tp (istep+1) 

z (istep+1) =z  (istep)  +  Qf*delt/ (cval*Tp (istep+1 )- (cvalinf* 
Tinfc) ) 

elseif (z (istep) . It . ze)  then 
z (istep+1) -ze 

Ql=subrhol*por*Aw*ze*subhsl 
if (imelt.gt.0)  Q1=0 . 

Tp (istep+1) =Tp  (istep)  +  (Qe-Qc-Ql) *delt/ (cval*ze) 
write  (10,*)  # ********  in  the  z  It  ze  spot z (istep) 

elseif (z (istep) .ge. (zt) )  then 

write  (10,*)  'z  is  out  of  here - ' 

Tp (istep+1) «Tp (istep) + (Qe-Qc) *delt/ (cval* (zt) ) 

z (istep+1) =zt 

imelt=0 

endif 

elseif (Tp (istep) .lt.Tmelt)  then 

cval= (piperho*pipecp*Ap)  +  (1 . -por) * (wickrho*wickcp*Aw)  + 
por*  <subrhos*subcps*Aw) 

Tp (istep+1) =Tp (istep)  +  Qe*delt/ (cval*ze) 

write  (10,*)  '  In  the  Tp  <  Tmelt  loop  and  cvalf=  ',cvalf 


elseif (Tp (istep) . It .Ttr . and.Tp (istep)  .gt.Tmelt)  then 
cval= (piperho*pipecp*Ap)  +  (1 . -por) * (wickrho*wickcp*Aw)  + 
+  por* (subrhol*subcpl*Aw) 

Tprop=Tp (istep) 
call  subprop 

Ql=subrhol*por*Aw*ze*subhsl 
if  (imeltm.eq.  1)  Q1=0  . 

Tp (istep+1) =Tp (istep)  +  (Qe-Ql) *delt/ (cval*ze) 
imelt^l 
imeltm=l 
endif 

istep=istep+l 
tval^delt* (istep-1) 
if (tval.gt .tmax)  goto  999 

write (11, 1030)  istep, tval, Tp (istep) , z (istep) ,hbare,hbarc 

goto  101 

continue 
format  (a20) 

format ('step  =  ',i4,f  Tp  =  ',f8.3,'  Toe=  ',f8.3,'  Tre=  ', 
f 8 . 3, '  T*r=  ' , f 8 . 3) 

format ('step  =  ',i4,'  Tp  =  ',f8.3,'  Toc=  ',f8.3,'  Trc=  ', 
f 8 . 3, '  T*c—  '  ,  f 8 . 3) 

format('Qs  =  ',el0.3,'  Qe  =  ',el0.3,'  Qc  =  ',el0.3,'  rho  = 
e7 . 2 , '  hfg  =  ' ,e7.2) 

format (lx,  i4 , f7 . 2, lx, f 8 . 3, lx, el2 . 6, lx, el 3 . 6, lx, el 3 . 6) 
end 


****************  Subroutine  airprop  ****************************** 
subroutine  airprop 

common /blkl/ao (8) , al (8) , a2  (8) , a3  (8 ) , a4  (8) , a5 (8) 
common/blk2/airk, airrhov, aircp, airmuv, airgam, airPr,  airR 
common/blk3/subrhov,  subrhol,  subrhos,  submuv,  subhfg,  subgam, 

+  subhsl,  subR,  subcpl,  subcps,  subcpv,  submul,  subsig,  subk 

common/blk4/piperho,  pipecp,  wickrho,  wickcp 
common /blk5 / Tpr op 

t*Tprop 


-  Using  curve  fit  data  - 

airrhov  =  2.45  -  5.72e-3*t  +  5 . 32e-6*t**2 .  -  1 . 70e-9*t**3 . 

aircp  *  0.9422267  +  1 . 926976e-4*t 

aircp  *  aircp*1000. 

airmuv  =  11.7858  +  0.6798*t  -3 . 8236e-4*t**2 .  +  1.1411e-7*t* 
airmuv  «  airmuv/ l.e7 

airk  «  -6.40514e-3  +  1.35212e-4*t  -  9 . 97069e-8*t**2 .  + 

+  3.73997e-ll*t**3. 

airPr  *  0.838509  -  6.98285e-4*t  +  9 . 82493e-7*t**2 .  - 

3.96097e-10*t**3. 
write (10,*)  t , airrhov 
end 


999 

1000 

1010 

+ 

1011 

+ 

1020 

+ 

1030 


+ 


c****************  Subroutine  pipprop 


****************** *******★★*★★ 


subroutine  pipprop 

common/blkl/ao (8) ,  al (8)  ,  a2  (8)  ,  a3 (8)  , a4  (8) , a5 (8) 
common/blk2/airk,  airrhov, aircpf  airmuv, airgam, airPr, airR 
common/blk3/subrhov, subrhol, subrhos, submuv, subhfg, subgam, 
+  subhsl,  subR,  subcpl,  subcps ,  subcpv,  submul,  subsig,  subk 

common/blk4 /piperho , pipecp , wickrho , wickcp 
common /blk5 /Tprop 

t*Tprop 


c - Using  curve  fit  data - 

piperho  *=  8328. 

pipecp  =  0.29418  +  6.51806e-4*t  -4 . 99916e-7*t**2 .  + 
+  1 . 63704e-10*t**3 . 

pipecp  =  pipecp*1000. 
c  pipek  = 

wickrho  =  piperho 


wickcp  =  pipecp 
end 

c****************  Subroutine  subprop  ****************************** 
subroutine  subprop 

common/blkl/ao  (8) , al  (8)  , a2  (8)  , a3  (8)  , a4  (8)  ,  a5 (8) 
common/blk2/airk,  airrhov,  aircp,  airmuv,  airgam,  airPr,  airR 
comraon/blk3/subrhov,  subrhol,  subrhos,  submuv,  subhfg,  subgam, 

+  subhsl,  subR,  subcpl,  subcps,  subcpv,  submul,  subsig,  subk 

common/blk4 /piperho, pipecp, wickrho, wickcp 
common /blk5 /Tprop 
common/blk6/iwf 

t*Tprop 


c  -  Using  B&K  Handbook  Least  Square  Data  - 

subPlog=ao  (1)  +  al(l)*t  +  a2(l)*t**2.  +  a3(l)*t**3. 
+  +  a4(l)*t**4.  +  a5 (1) *t**5 . 

Press*exp (subPlog) 

subrhol=ao (2)  +  al(2)*t  +  a2(2)*t**2.  +  a3(2)*t**3. 
+  +  a4(2)*t**4.  +  a5  (2)  *t**5  . 

c  **  B&K  gives  nu,  not  mu  ....  (in  the  table  of  values) 

submul=ao(3)  +  al (3) *t  +  a2(3)*t**2.  +  a3(3)*t**3. 
+  +  a4(3)*t**4.  +  a5 (3) *t**5 . 

subrhov«ao (4)  +  al(4)*t  +  a2(4)*t**2.  +  a3  (4) *t**3 . 
+  +  a4(4)*t**4.  +  a5  (4) *t**5 . 

subrhov=exp (subrhov) 

c  **  B&K  gives  nu,  not  mu  ....  (in  the  table  of  values) 

submuv=ao(5)  +  al (5) *t  +  a2 (5) *t**2 .  +  a3(5)*t**3. 
+  +  a4(5)*t**4.  +  a5  (5) *t**5 . 

subsig=ao{6)  +  al(6)*t  +  a2(6)*t**2.  +  a3  (6)  *t**3 . 
+  +  a4(6)*t**4.  +  a5 (6) *t**5 . 
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subhfg=ao (7)  +  al  (7) *t  +  a2(7)*t**2.  +  a3(7)*t**3. 
+  +  a4  (7) *t**4 .  +  a5 (7) *t**5 . 

subk-ao  (8)  +  al(8)*t  +  a2(8)*t**2.  +  a3  (8) *t**3 . 
+  +  a4 (8) *t**4 .  +  a5  (8) *t**5 . 

if(iwf.eq.l)  then 


c - -  Using  curve  fit  data  for  Sodium - -  — 

subrhos=-0 .22127*t  +  1033.45338 

subcps  =  0.5959051  +  4.838459e-3  *t  -  1.450932e-5  *t**2. 
+  +  1 . 73901e-8  *t**3. 

subcps  «=  subcps  *10  00 . 

c - Need  to  find  actual  data - 

subhsl=113044. 

c  units  of  above  are  J/kg 

subgam=l . 67 


c  specific  heat  data  from  Incropera  and  DeWitt,  Table  A. 7 
subcpl^l . 3 
subcpl=subcpl* 10 00 . 

write  (9, 1010)  t, press , subrhol, submul, subrhov, submuv, subsig, 
+subhfg,  subk, subrhos, subcps, subgam, subcpl 
endif 

if ( iwf . eq . 2 )  then 

c - Using  curve  fit  data  for  Potassium - 


subrhos=-0 . 1 6743*t  +  914.46661 

subcps  =  0.5073591  +  1.357728e-3  *t  -  2.079707e-6  *t**2. 
+  +  9 . 896415e-10  *t**3. 

subcps  =  subcps* 1 000 . 

c - Need  to  find  actual  data - 

subhsl=61127 .28 

c  units  of  above  are  J/kg 

subgam=l . 4 


c  specific  heat  data  from  Incropera  and  DeWitt,  Table  A. 7 
subcpl=0 .75 
subcpl^subcpl^OOO . 

write (9,1010)  t, press, subrhol, submul, subrhov, submuv, subsig, 
+subhfg, subk, subrhos, subcps, subgam, subcpl 
endif 

1010  format  (f  8 . 3, 10  (lx,  e9 . 2)  ,  lx,  f3 . 1,  lx,  f7 . 2) 
end 
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